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\ A mesoscopic or coarse-grained approach is presented to study thermo- 

Q>^ I capillary induced flows. An order parameter representation of a two-phase 

O ' binary fluid is used in which the interfacial region separating the phases natu- 

rally occupies a transition zone of small width. The order parameter satisfies 
the Cahn-Hilliard equation with advective transport. A modified Navier- 
Stokes equation that incorporates an explicit coupling to the order parameter 
d ' field governs fluid flow. It reduces, in the limit of an infinitely thin interface, 

to the Navier-Stokes equation within the bulk phases and to two interfacial 
forces: a normal capillary force proportional to the surface tension and the 
mean curvature of the surface, and a tangential force proportional to the tan- 
gential derivative of the surface tension. The method is illustrated in two 
cases: thermo-capillary migration of drops and phase separation via spinodal 
decomposition, both in an externally imposed temperature gradient. 
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I. INTRODUCTION 



The study of multi-phase flows leads to a classical moving boundary problem in which the 
equations governing fluid motion are solved in each phase, subject to boundary conditions 
specified on the moving boundaries. In the classical case, the boundary of separation is 
assumed to be an idealized boundary without any structure. For viscous flows, the velocity 
field is continuous across the boundary, whereas the normal component of the stress tensor 
is discontinuous if capillary forces are allowed for. Implicit in this formulation are, of course, 
assumptions that local thermodynamic equilibrium obtains and that the typical length scales 
of spatial structures and flow are much larger than the scale of the physical boundary 
separating the phases. In most cases of interest to fluid mechanics the scale of the flow is set 
externally by macroscopic stresses on the system, whereas the length scale of the interface 
is of microscopic size and set by the range of the intermolecular forces in the fluid. Under 
these conditions, the transition region between the two phases can be approximated for all 
practical purposes by an ideal, mathematical surface of discontinuity. 

A notable exception to this rule concerns flow in a fluid near or at its critical point. 
There the surface of separation between both phases becomes arbitrarily diffuse, and the 
flow within the boundary itself needs to be explicitly modeled and resolved. The so-called 
Model H (following the nomenclature of Hohenberg and Halperin |I|]) was introduced to 
study order parameter and momentum density fluctuations near the critical point, as well 
as their interaction, at a coarse-grained or mesoscopic scale p|. Additional (reversible or 
Hamiltonian) terms were added on phenomenological grounds to the Navier-Stokes equation 
and to the equation governing the relaxation of fluctuations in the order parameter appro- 
priate for the phase transition under consideration. Later studies based on the same model 
have addressed critical fluctuations of simple fluids and polymers |^,|| under external 
shear, and more recently (and more closely related to our study) hydrodynamic effects on 
spinodal decomposition in a binary fluid 0J|] . In spinodal decomposition, a binary system is 
quenched from a homogeneous disordered state to a region of the phase diagram in which the 
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homogeneous state is linearly unstable to fluctuations in a wide band of length scales. (For a 
review, see, e.g., 0.) Some further details on the general principles behind a mesoscopic de- 
scription of equilibrium properties of bulk phases and of interfacial structures can be found 
in readily available texts and monographs (see, e.g., Refs pO|-p^.) Implications of such 
a level of description on the equations governing fluid flow are discussed in the appendix, 
where the modifications of the non-dissipative part of the stress tensor embodied in Model 
H are obtained. This "Hamiltonian" or reactive part of the stress tensor does not contribute 
to the entropy increase. Additional treatments based on irreversible thermodynamics and 
continuum mechanics can be found in Refs. |jl5|,|l6|. 

The usefulness of a coarse-grained approach away from a critical point is less well es- 
tablished. Previous work along these lines includes Ginzburg-Landau models used in the 
study of the kinetics of phase separation (e.g., the Cahn-Hilliard equation) 0, or phase 
field models of solidification In both mathematical surface of separation 

between the two phases becomes a transition region of small but finite width, across which 
all magnitudes change continuously. The essential ingredient in both cases is the assumption 
that the relevant thermodynamic potential density (entropy in closed systems, Helmholtz 
free energy if the system is held at constant temperature) is a function not only of the usual 
thermodynamic variables and of the order parameter, but also of the spatial gradient of the 
latter. In the case of a binary fluid, for example, the entropy density is assumed to depend 
on both the mass density of solute and its spatial derivative. As a consequence, the chemical 
potential of the solute is no longer given by the derivative of the local entropy with respect 
to solute density. 

From a theoretical point of view, it is plausible to assume that modes describing the 
motion of boundaries separating bulk phases, even away from a critical point, are the long- 
lived modes of the dynamics of the system, and hence gradients of the order parameter can 
be considered as additional arguments upon which the thermodynamic potential depends, 
at least on the time scale over which the motion of the boundaries is observable macro- 
scopically. Furthermore, one would also like to satisfy the consistency requirement that the 
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coarse-grained model reduces to the classical macroscopic description in the limit of a small 
transition region (compared to the scale of relevant structures). This requires, in general, 
relating macroscopically measurable quantities (e.g., the surface tension of the boundary) 



to phenomenological parameters that enter the mesoscopic model ||20|| . 

One of the chief advantages of a mesoscopic approach is that it allows the study of fluid 
phenomena that lie outside the classical continuum formulation, for example, the break 
up and coalescence of fluid domains. More generally, by incorporating explicitly into the 
model dissipation at short length scales (on the order of the interfacial thickness or order 
parameter correlation length), the model allows the study of other physical situations that 
may be influenced by phenomena at that scale. Other examples that may fall in this class 
include motion of contact lines, flow near solid walls and slip at a microscopic or mesoscopic 
scale. 

From a purely computational point of view the model used can be viewed as an extension 
of a class of numerical methods often used to study the classical moving boundary prob- 
lem, namely, those that phenomenologically introduce additional dissipation at short length 
scales. At least in principle, however, the method that we describe in this paper differs from 
those approaches (e.g., methods based on the introduction of ad-hoc "artificial dissipation" 
or "artificial viscosity"), since, in the present case, short length scale dissipation is intrin- 
sically part of the model, and is related to the order parameter variable that describes the 
phases. Also from a computational point of view, the approach that we follow is similar to 
the so-called immersed boundary or diffuse interface methods introduced for the study of 
multi-phase flows [pT| -|25[| . The order parameter which we define below plays a role analogous 



to the "color function." However, it is not an auxiliary field introduced for computational 
convenience, but rather a thermodynamic variable in its own right, controlled by a local free 
energy density. As a consequence, for example, the surface tension of this model cannot be 
fixed independently, but is completely determined by the choice of thermodynamic poten- 
tial and the parameters contained therein. This will reveal some interesting features to be 
expected at high thermal gradients. 



In Section || we introduce the coarse-grained model used in our study. We concentrate 
in this work on thermo-capillary flows and hence exphcitly discuss how to couple the equa- 
tion governing fluid motion and relaxation of the order parameter. In this paper we will 
consider the case of a constant, externally imposed temperature gradient, although it is 
straightforward to extend our calculations to incorporate a fluctuating temperature fleld. 
We also discuss in this section the limit of a thin interface and the boundary conditions 
that result from our model in that limit. As illustrations, we apply the numerical approach 



to two problems of interest. In Section |IT| thermo-capillary motion of small drops under 
large thermal gradients is considered, while in Section |3 we study phase separation via 
spinodal decomposition in the presence of a uniform temperature gradient. Section is 
reserved for concluding remarks. An appendix provides additional details on the derivation 
of modiflcations to the Navier-Stokes equation introduced through the coarse-grained level 
of description. 



II. COARSE-GRAINED MODEL 
A. Formulation 

Consider an incompressible and Newtonian binary fluid at a temperature below its phase 
separation critical point. Let ip be the order parameter density appropriate for the un- 
mixing transition, chosen to be = in the disordered phase above the transition point, 
and symmetric and equal to -^^Peq below, [ip can be thought of as c — Cc, where c is the mass 
fraction of one of the species and Cc its value at the phase separation critical point.] One 
further makes the reasonable assumptions for binary liquid systems that the shear viscosities 
of both phases are equal, and that the dependence of the mass density on the order parameter 
is weak and can be neglected. In effect, we consider here a "symmetric" model in which all 
bulk properties of both phases are equal. These restrictions, which conveniently dramatize 
interfacial phenomena in the cases which we study here, can be easily removed. 
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If V is the velocity of an element of fluid, conservation of mass and momentum lead to, 

V-w = 0, (1) 



= -Vp + r^V^tT+ZiVy?, (2) 

where /i = 6J-'/6(p is the chemical potential conjugate to the order parameter (f {S/6(f stands 
for functional or variational differentiation), p and rj are the density and shear viscosity of 
the fluid respectively, and p is the pressure field. The last term in Eq. is non-standard 
in macroscopic approaches and incorporates capillary effects as discussed further below and 
in the appendix. As a reference, we note that the Navier-Stokes equation is also modified 
in Immersed Boundary Methods by adding a term of the form, 

(^ann + 6 {x - xs) , (3) 

where is the instantaneous location of the boundary. Here a is the interfacial tension 
(independently prescribed), n and s are the unit vectors along the normal and tangential 
directions respectively, and n is the mean curvature of the interface. As we further discuss 
in Section |lllj| , the term fiVif in Eq. in the limit of a thin interface reduces to the term 
used in Immersed Boundary Methods, but with a now determined by the coarse-grained 



free energy functional JF (see, e.g., Refs |1T0|-|14|) that, for simplicity, is chosen to be of the 
form, 

J dV [^\V^\' + (4) 

with, 

= -^^^ + j^'- (5) 

The integration extends over the entire system (both bulk phases and interfacial regions), and 
K, r, and A are three phenomenological coefficients as yet unspecified, other than requiring 
K,X > 0. For the model defined by Eqs. (§) and (^ the chemical potential is given by 



II = —KV'^ip — TLf + Xif^. As is well known, this free energy qualitatively describes a single 
homogeneous phase over a range of parameters (r < 0) yielding to two-phase coexistence 
in another region (r > 0). The parameter r is therefore a measure of distance below the 
phase separation critical temperature and is taken to be proportional to that temperature 
difference, while the parameter A, at this level of approximation can be related to the inter- 
particle potentials (or to a virial coefficient). As noted below, the parameters A and K can 
be determined from equilibrium measurements. The order parameter is further assumed to 
satisfy a modified Cahn-Hilliard equation, 

^ + V.{^v) = MVV (6) 

to allow for advective transport of if. M is a phenomenological mobility coefficient of 
microscopic origin, which, in a real system, could be inferred from mutual diffusion and 
order-parameter susceptibility measurements away from the critical region. 

Equations (|l])-(0) and (|)-(|) completely specify the model. We have imposed the fol- 
lowing boundary conditions at the edges of the fluid domain: {7 = 0, = ipeg{T), and 
dip/dn = dipeq(T) / dn, i.e., the value of the order parameter near the boundary is deter- 
mined by the imposed temperature field there. Note that with this choice of boundary 
conditions, / dV(f is not a strictly conserved quantity. The boundary conditions required 
to enforce strict conservation of order parameter are dip/dn = and d'V'^ip/dn = 0. In 
this case small boundary layers of (f (of thickness of the order of the correlation length ^) 
would develop because ip in the bulk (and away from interfaces) is approximately equal to 
iPeg(T), and therefore would not satisfy this latter set of boundary conditions. Our choice 
of boundary conditions is motivated by computational simplicity so that we do not have to 
resolve additional boundary layers at the edges of the domains. In any event, the spatial 
integral of changes very little during the course of the numerical solution, and hence this 
choice has little effect on the dynamics in the bulk phases. 

From Eq. (^ one finds ipeg = with ±ipeq as the two coexisting solutions. Equation 

(P) also admits a nonuniform, "kink" stationary solution for -u = and constant r > 0, 
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ipoiu) = ipeqtOJoh-j^, (7) 

with boundary conditions (p — * iy^gg as m ^ ±00. The coordinate u is along an arbitrary 



direction and ^ is the mean field correlation length, ^ = \jK/r. This stationary solution 
represents the coexistence of two equilibrium phases with ip = ±(peq, across a transition 
region located at m = and of width C,- In what follows, we will refer to the regions 
\u/^\ ^ 1 as the bulk phases, and \u/C,\ ~ 1 as the interfacial or transition region. In the 
region of parameters of Eqs. (|) and (|]) that we explore, the surface area occupied by the 
bulk phases is much larger than that of the interfaces. 

The excess free energy density associated with Eq. compared to the free energy of 
either uniform phase is given by |T^, 

We introduce dimensionless variables by choosing ^ as the scale of length, K/Mr"^ as 
the scale of time, and ipeq as the order parameter scale. In dimensionless units, the model 
equations reduce to, 

V-^T=0, (9) 
Re^ = - Vp + V^v + C (- V V -V + V"^) V</^, (10) 



and. 



^ + V ■ {p>v) = (- W -^ + ^'). (11) 



Two dimensionless groups remain; one of them, Re= Mr/u {u = rj/p is the kinematic 
viscosity) gives the ratio between order parameter diffusion and momentum diffusion due to 
viscosity. In the calculations that we describe in this paper, the system is in the overdamped 
limit of Re~ 0. The other, C= 3a^/2r]Mr, plays the role of a capillary number. Flows in 
our study are entirely driven by surface tension. In the overdamped limit the velocity scale 
of such flows is V ^ a/rj where we have assumed that the scale of the flow is the same as 
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the scale of the typical domains of the two phases. Therefore, the characteristic time for an 
element of fluid to be advected a distance equal to the interfacial thickness is = ^rj/a. 
The diffusive time scale given above is = ^'^ /Mr; hence C oc t^/tc,. 

In order to couple the previous equations to a slowly varying temperature field, we make, 
as noted above, the traditional identification of the parameter r in Eq. (|^) as being a linear 
function of temperature r{x) oc Tq — T(x), where To is a constant playing the role of the 
consolute or critical temperature. In the case that we study, the temperature field is fixed 
and remains constant throughout the calculation. If one wishes to incorporate a fluctuating 
temperature field, the model equations have to be supplemented with the equation of energy 



conservation [Q. In dimensionless units, we consider, 

^ = -VV - + <^', (12) 

where the dimensionless local temperature variable r(y) = ro + ay, with a the dimensionless 
temperature gradient. We identify the direction as the direction in which the imposed 
temperature varies. For our coarse grained analysis to make sense we require that temper- 
ature variations are small over small distances (of the order of the interfacial thickness), 
namely, in dimensionless form, |a| <^ 1, but perhaps large over distances of the order of the 
system size, and possibly even over a domain of one of the phases. 



B. Sharp interface limit 

We next present a heuristic analysis to illustrate that the coarse-grained model introduced 
reduces to the classical continuum description in the limit in which the width of the interfacial 
region is much smaller that the size of the bulk domains. We introduce a local orthogonal 
system (s, u) such that a given point in space is given by r = R{s) + uh, where R{s) 
describes the location of the level set ^{r) = (s is the arclength along it), n is the local 
normal pointing towards the (+) phase, and u is the distance along the normal direction. 

In the bulk regions (Im/^I ^ 1) yuVy^ is of order higher than linear in the gradients and 
therefore negligible. One then recovers the standard Navier-Stokes equation. The Cahn- 
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Hilliard equation can be likewise linearized around ipeq in the bulk regions, yielding the 
standard diffusion equation with advective transport. 

In the interfacial regions ~ 1)? the term /iVy? does become large. In the limit 

of gently curved interfaces, ^ 1, where k is the mean curvature of the interfacial 
region, and for small temperature gradients, |a| -C 1, the field relaxes on a fast time 
scale to essentially the stationary solution ipoiu) as in Eq. (|^ |^. One can develop a 
multiscale expansion of the field ip = ip{S, u, U) with S, U the slowly varying coordinates in 
the tangential and normal direction as set by the slow variation of the temperature field, and 
u the fast variation in the normal direction as determined by the free energy functional. We 
further assume - and this remains to be proven more rigorously - that ip relaxes quickly to its 
local equilibrium profile given by the local value of the temperature, but there remain slow 
gradients of (f, since this approximate solution is no longer a solution of 5T /5ip = constant, 
throughout the system. The surface tension in Eq. (P) can then be written as, 

/ + 00 
duAF{u,U,S), (13) 
-oo 

where F = ^K\'V(f\'^ + f{(f) and AF is the free energy (density) difference between the mean 
field value of F at the local value of in a configuration with a two-phase interface and the 
free energy of the bulk phases ipeq, at the local temperature. Under the assumptions stated 
above, this expression can be used to define a slowly varying surface tension (j{S) along the 
interface, since the magnitude of the order parameter will be the local equilibrium value 
at the local temperature. Then, 




where we have used the fact that the metric tensor element g22 = ^ ' ^ = (1 + uk)"^ ^ 
1 in the limit of gently curved interfaces considered. (Note that k < for a sphere). 
The chemical potential appearing is appropriate for an interfacial configuration and is non- 
vanishing for a gently curved (or indeed a flat) interface in a temperature gradient. Hence 
the tangential component of fiV(f gives rise to a tangential surface force which equals the 
tangential derivative of the surface tension. 
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We can likewise study the normal component by recalling that, 



dufi 



dip 
du 



du 



-oo 
oo 

~ ' du 

' — oo 



du 



2 "^^^9^ 



du 
dip d'^p 



^ dui^^ 
-oo I du 



dS^ 

^ ^ d'^ip dp 
oo dS'^ du 



dp 
du 



(15) 



The first term gives again the surface tension a, and the second term would remain in the 
sharp interface limit due to the slow variation of p with temperature. 

Therefore in the case of a thin interfacial region, the additional term in the Navier-Stokes 
equation is equivalent - under the stated assumptions of a slowly varying temperature field 
and gently curved interfaces - to a tangential surface force proportional to the tangential 
derivative of the surface tension, and to a normal force proportional to the mean curvature 
of the interface. 



C. Numerical method 

We next discuss the numerical algorithm used in the actual computations. We restrict 
our attention in this paper to two dimensional flows. Since the velocity field is solenoidal, 
it is convenient to introduce the stream function, 

-Vx(#.) = |S-f3. (16) 

where the velocity field is defined in the (x, y) plane, and k is the unit normal perpendicular 
to that plane. By taking the curl of Eq. ( [10|) with Re = we find. 



vV + c [V [y'p) X v^J ■ /c = 0. (17) 

The flow field can be found by solving a biharmonic equation for the stream function, subject 
to the boundary conditions that both the stream function ip and its normal derivative dip/dn 
vanish on the external boundaries of the system. 
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Concerning the equation for the order parameter, Eq. ([TT|), we follow Ref. jS^] and use a 
backward implicit method which is first order in time, 

"^^^ ^^At~ '^^^^ + ^ ■ ^"^^^^^^^ + ^^^^ = + ^''^'^^ ^^^^ 

with £ = V^H- V^. The convective term has been kept in conservative form, and the velocity 
field retained in the discretization is at time t. The boundary conditions that we use for the 
order parameter are (p = (feq{T) and dip/dn = d(peq{T)/dn on the external boundaries. A 
Gauss-Seidel iteration is used to solve Eq. (^). First, one considers an "outer" iteration, 

ip{t + At) c:^ Pk+i = <Pk + Sk with ipo = ip{t). (19) 

Substituting ip(t + At) in Eq. ( pTSD and linearizing in the outer correction field 5k{x, y) yields 
an equation for the outer correction with a residual that is a function of (pk- Successive 
iterates converge to a solution when both the residual and go to zero simultaneously. 
The outer iteration is performed simultaneously on Eq. ( P^ ) and (|T^. In practice, an 
"inner" iteration is set up to solve for 5^, by assuming, 

4 - 5k,m+l = Sk,m + Vm, (20) 

where rjm is the inner correction. The variable coefficient biharmonic operator acting on 77^ 
is then approximated by a constant coefficient biharmonic operator. 



where a = 1 — 3 < pi > and b=l/At — 6< pk^'^Pk >, and < • > stands for spatial averages 
over the entire system. Then both Eq. (p^ and (|2TD are solved with a fast biharmonic solver 



see ref. [30| or [31] for additional details). 



III. THERMO-CAPILLARY MOTION OF DROPS 



In this section we briefiy summarize application of the computational method introduced 
in Section |I| to the thermo-capillary fiow of a drop of one phase in the background of its 
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coexisting partner. As we have noted above, the main strengths of the coarse grained 
technique are that (i) it allows natural tracking of the interface separating the two phases, 
(ii) it naturally encompasses topology changes, such as the merging of two droplets, and (iii) 
it can include phenomena from the scale of the thermal correlation length up to macroscopic 
while providing a qualitatively correct picture of behavior near criticality. Accordingly, to 
separate effects from those which are normally treated at the macroscopic level, we turn to 
relatively small drops and relatively large temperature gradients. 

To illustrate the method we consider the thermo-capillary flow of a single drop. Initially 
the radius of the small drop is fixed at i? ~ 10^, where ^ is the thermal correlation length at 
our reference temperature tq = 1, which is typically the value of the temperature parameter 
at the cold end. This size drop is about as small as we can go within the coarse grained 
description and retain some expectation that the qualitative results will translate to the 
macroscopic scale. However on a macroscopic scale, we consider temperature gradients which 
are large, so that the dimensionless number is R\WT\/T ~ 0.01. Hence, over a correlation 
length, the temperature change is typically 0.1%, allowing us to retain the coarse grained 
approach, but over a drop, the change is on the order of 1%. By macroscopic standards 
this is large. Consider, for example, a nucleating 1/im drop in a binary liquid system in a 
temperature gradient of 1 C/mm at roughly room temperature. The dimensionless ratio is 
three orders of magnitude smaller. 

Our first diagnostic is the drift velocity of such a droplet as a function of the temperature 
gradient. In all figures to follow, the high temperature side is at the bottom. First a sample 
plot of the droplet "center of order parameter" (defined by analogy with the center of mass) 
as a function of time is shown in Fig. |l|, indicating motion with constant velocity. We have 
repeated the calculation for a number of values of the temperature gradient a, and the 
results are shown in Fig. 0. At low values of the gradient, we observe linearity ||32|,Q of 



velocity, but at our higher values there is a reproducible reduction as seen in the figure. 
Furthermore, at the higher values of the gradient the velocity of the drop is dependent on 
the temperature itself, as well as the gradient. These features are associated with the fact 
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that this model, which can (quahtatively) apply near the phase separation critical point 
as well, has an equilibrium order parameter depending on the dimensionless temperature 
parameter as r^/^ ~ (Tq — T)^/^, as well as a nonlinear dependence of the interfacial tension, 
(J r^/'^ (To — T)'^/^. These effects, for fixed gradient, become relatively more apparent 
near the high temperature side. 

Despite the fact that the drop velocity for high gradients reflects nonlinearities inherent 
in the model, we have found that, for small enough temperature gradient, the velocity scales 
quite linearly with drop radius over a limited range (a factor of three) accessible in these 
preliminary studies. For larger gradients, the explicit dependence of both order parameter 
and surface tension on temperature eventually affects the proportionality as seen in Fig. ^ 

Even at the relatively high gradients (by macroscopic standards) used we observe little 
distortion of the droplets. Discussion of this issue, further details on the flow fields and 
analytic analysis of the droplet migration analysis based on the coarse grained equations is 
beyond the present scope. 

We have finally done a qualitative study of drop coalescence in a temperature gradient. 
A sequence of configurations illustrating the motion of two drops is shown in Fig. ^ Co- 
alescence involves a topology change, which is naturally handled within the coarse-grained 
approach. An analysis of the kinetics of droplet coalescence, and droplet detachment and 
attachment to boundaries, is a potentially rich area but also beyond the present scope. 

IV. SPINODAL DECOMPOSITION IN A TEMPERATURE GRADIENT 

In order to study whether the method can describe more complex flows, we have studied 
spinodal decomposition of a binary fluid, in two spatial dimensions, and under an imposed 
constant temperature gradient. In dimensionless units, the size of the rectangular fluid 
domain studied is a = 800 and b = 200 (in the gradient (y) direction). We have used an 
evenly spaced grid with m = 1024 nodes along the x direction, and n = 256 along the y 
direction. The results presented involve a dimensionless temperature gradient a = —0.004 
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along the y direction, and include up to t = 100. The initial condition for the order parameter 
at each point is a set of uniformly distributed random numbers in [—0.01,0.01]. The time 
step used in the numerical integration is dt = 0.25. 

Figure ^ shows an example of evolution of the system for three different instants of time: 
t = 6.35, 50 and 100. The order parameter (between -1 and 1 in the dimensionless units 
we are using) is shown in grey scale. The characteristic spinodal pattern emerges, with an 
intensity and domain size that is a function of the local temperature of the system. The 
temperature parameter r as introduced in Eq. @ is known to couple extremely weakly to 
the phase separation process. It is not expected to change the algebraic form of the growth 
law for the domains (i.e., the power-law growth in time of the characteristic domain size), 
but rather it may change the amplitude. Two questions naturally arise for a slowly varying 
temperature, namely, whether the temperature gradient introduces any anisotropy in the 
characteristic length scale of the pattern, and whether an effective growth law amplitude 
exists that slowly changes in space, as the temperature of the system changes. Roughly 
speaking, the question is whether in the presence of a slowly varying temperature, phase 
separation proceeds temperature "slice" by "slice" as in ordinary spinodal decomposition at 
that final temperature. Directional anisotropy appears to arise in the purely diffusive Model 
B (conserved order parameter without hydrodynamic coupling, fl]]) at early to intermediate 
times (in the form of a lag in the growth parallel to the gradient), and it carries over into the 
late stage growth stage. |35| Preliminary results over a small ensemble of independent initial 
conditions (five) indicate that, contrary to the situation in Model B, there is no appreciable 
anisotropy in the characteristic length scale of the domains. It is interesting to note that 
even if thermo-migration effects themselves are small for dimensionless times t < 100, the 
hydrodynamic coupling and complex flows appear to wash out any prominent anisotropy. 
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V. CONCLUSIONS 



We have implemented a numerical algorithm to solve the coupled Cahn-Hilliard and 
extended Navier-Stokes equations (Model H) in a temperature gradient. The backward 
implicit method is unconditionally stable, and with moderate computing effort has allowed 
us to study reasonably large system sizes for a reasonable amount of time. Two particular 
examples of thermo-capillary induced flow in two dimensions have been studied: the motion 
of droplets of one phase in the background of its coexisting partner in a temperature gradient, 
and phase separation via spinodal decomposition, also in a temperature gradient. In two 
dimensions a stream function representation has been used so that an efficient biharmonic 
solver handles both the order parameter and flow equations together. In the first case, 
the results given by this method agree with classical macroscopic calculations in that the 
migration velocity of the drop is proportional to both the imposed temperature gradient 
and the radius. For large gradients, the velocity depends nonlinearly on the temperature 
gradient, as it is to be expected given the dependence of the miscibility gap and surface 
tension on temperature in the modeling of the coarse grained free energy. 

The chief advantage of the method is that it allows the study of the fiow on length 
scales which are not accessible to classical macroscopic methods. One such case involves 
the coalescence of drops and the concomitant topology change of the interfaces. We have 
presented qualitative results involving the coalescence of two drops induced by thermo- 
capillary migration. Of course, the mesoscopic or coarse-grained model introduced specifies 
a dependence of the fiuid variables at such length scales, and this requires the introduction 
of a phenomenological free energy at that coarse-grained scale. This coarse-grained free 
energy is a function of a few phenomenological parameters (imagined to depend smoothly 
on the temperature and other experimentally accessible variables) and couples through the 
associated chemical potential to the velocity field. The requirement that the model equations 
lead to the classical, macroscopic description in the limit of a small interfacial region allows 
one, in principle, to determine some of the phenomenological parameters (at least within a 
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mean field context). 

We have also argued that the model introduced reduces to the classical boundary condi- 
tions at the two-phase interface: namely, a normal force given by the surface tension of the 
interface and its mean curvature, and a tangential force equal to the tangential derivative 
of the surface tension. These boundary conditions follow in the limit in which the imposed 
temperature field does not change appreciably over the scale of the interface. For these 
arguments we have used the plausible assumption (consistent with observations from simu- 
lations) that the order parameter within the two phases and in the interfacial region relaxes 
quickly to its local equilibrium value determined by the local curvature of the interface and 
the local temperature. A more systematic analysis to elucidate this point is clearly needed. 
This analysis can also lead to a systematic cataloging of corrections to the macroscopic 
boundary conditions. 

An additional advantage of the method is that it allows the study (at least quahtatively) 
of very small drops and large gradients. In fact, this method is not likely to be competitive 
with classical methods in situations in which gradients are weak and drops are large, as 
compared to the scale of the interface. We remark, however, that classical methods cannot 
handle topology changes that are controlled by physical processes on short or intermediate 
length scales, such as the order-parameter correlation length. Finally, there are several 
classes of problems in which the behavior at short lengths scales needs to be resolved since 
it eventually determines the fiow at macroscopic scales. Examples include the motion of 
contact lines, discontinuous velocity fields at boundaries, and slip and thermal fiuctuations 
at the mesoscopic scale. 
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APPENDIX A: COARSE-GRAINED FLUID MECHANICS 

In this Appendix, a derivation is presented of the modifications of the Navier-Stokes 
equation that are introduced if the description of the system goes beyond local thermody- 
namics. In other words, if the density or order parameter varies on short length scales so that 
a square gradient or similar description is appropriate (see, e.g., Refs [p^OHI^), the Navier- 
Stokes equation contains new terms in the dissipationless (reactive or Hamiltonian) part of 
the stress tensor. To find these terms we can consider an ideal fluid. For this derivation we 
use a Hamiltonian (canonical) formalism since, we believe, it is simplest. A more detailed 



discussion of the Hamiltonian formalism is contained in the book by Zakharov et al. |^ 
which we draw on; a discussion of variational principles for an Euler fluid is contained in 
the monograph by Serrin. B? 



Following Zakharov et al. we consider an inviscid fluid with 

dtp + V- (pv) = 

^tX^ + v■Vxi = (Al) 

where {/(r, t) is the Eulerian velocity field, p(r, t) is the mass density and the scalar functions 
Xi, i = 1,2,..., represent any advected quantities carried with the flow and satisfying 
dx/dt = 0. The entropy per unit mass satisfies ds/dt = in isentropic flow so that we may 
treat s = constant. This allows the pressure to be expressed as p = p{p) and the local part 
of the energy per unit volume may be taken as e(p, s) —>■ e(p). For a single component fluid 
the function x represents an advected velocity circulation. Later when we consider a binary 
fluid, in the absence of diffusion, the order parameter ip, representing the mass fraction of 
one species, will be such a quantity. As usual the local part of the enthalpy per unit mass 
is h = de/dp. The energy of the fluid is taken to be 



n = J llpv' + e{p)]df. (A2) 



Following Ref. I^G] we show that standard methods reproduce the Euler equation. Then we 
modify Eq. ( |A2| ) to go beyond local thermodynamics to find the new contributions to the 
stress tensor. 

A Lagrangian is defined in a standard fashion introducing Lagrange (undetermined) 
multipliers for the conservation of p and 

C = J {^^pv'-eip) + ^[dtp + V-ipv)] (A3) 
-X[da + {v-Vx)]}dr (A4) 

The total action S = J Cdt is presumed to be an extremum with respect to fiows v which 
yields 

^=V$ + -Vx (A5) 
P 

Variation with respect to the multipliers $ and A recovers the conservation conditions in 
Eqs. ([Al|). Substitution of the velocity back into the Lagrangian yields the hamiltonian 
(canonical) description 

C = J [^dtp-Xdtx]df-n (A6) 

where Ti. is expressed in terms of the conjugate pairs p, $ and A, x- The Hamiltonian 
equations of motion are 

dtp = 6n/6^ 

= -6H/6p (A7) 
dtX = 6n/6x 
dtx = ~6H/6X 

It is straightforward to show that the second and third equations yield 
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dt<l> + V ■ V<l> - ^v^ + h = (A8) 
dtX + V- (Xv) = (A9) 

Applying V to the first yields the familiar Euler equation for the velocity change 

dtv + v-Vv = -p'^Wp (AlO) 

using the local thermodynamics dh/dp = p~^. This is easily generalized to include external 
potentials yielding additional body forces. 

Now if the physics dictate that one must go beyond local thermodynamics with, for 
example, a square gradient, the Hamiltonian ( |^ ) is modified to read 



n 



lpv'+<p)+lKip)ivpr 



dr. (All) 



Following all the same steps yields the modified Euler equation 



dtV + V ■ Vv = p ^ 



-Vp + pV{V-KVp)-^K'{p){Vpf 



(A12) 



Terms other than can be thought of as additional body forces owing to short length 
scale variations of the density. The dissipationless part of the stress tensor is appropriately 
modified. In the most general case, any functional g{p, Vp) can be added to the local energy 
in Eq. (|A11|) , and the stress tensor is modified accordingly. The modification of the stress 



tensor for a single-component fluid at the square-gradient level was derived by Felderhof ||38 
using a Lagrangian (particle) description of the mechanical equations. 

Further discussion of the generality of variational formulations of inviscid hydrodynam- 
ics are contained in the Refs. [^0. The addition of one additional advected quantity, 
representing the "conservation of the identity of fluid particles" will allow every flow to be 
represented as an extremal for the (Herivel-Lin) variational principle. ||37 



Now we turn to the more interesting case of a two-component fluid. Background is con- 
tained in the text of Landau and Lifschitz |^ . The total momentum flux is pv, where p is 



the total mass density. The density of one of the species is pc where c is the (dimensionless) 
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mass fraction, which plays the role of the order parameter (denoted (p in the body). Now the 
local energy per unit volume becomes e(p, c, s), and for isentropic flow of an ideal fluid we 
may neglect the specific entropy. In the absence of diffusion we have dtc + v ■ Vc = so that 
we may consider c as an advected quantity, denoted generally by x above. For simplicity 
(and realistically for typical binary liquid systems) we consider modifications of local ther- 
modynamics due to short length-scale concentration variations only. In other words, for the 
present purposes we assume, as is normally the case, that any density variations are on scales 
large compared to the order-parameter correlation length, ~ ^, so that a square-gradient 
contribution for the total density is not required. This assumption may be generalized in 
special circumstances. 

Hence, the total energy of the system is taken to be 



dr. (A13) 



In Hamilton's equations (A7) the mass fraction c now plays the role of the advected x as 
noted. Now STi/dc contains new terms, and the A equation thus becomes 

dt{-X) = V ■ (At/) + V ■ KS/c - ]^K'{c){Vcf - dt/dc . (A14) 

Note that dTi/dp remains the same since we have not included Vp as noted above. From 
Eq. (^) one evaluates the acceleration as 



dtv = Vdt^ + {\/p)VdtC + {dt{\/p)) Vc (A15) 
which yields, using Eqs. (^) and (|A14|) , 

^ Vc. (A16) 



dtV + {v ■ V)v = —p ^\/p + p ^ 



de/dc - (V ■ KVc) + ^K'{c)(ycf 



We have also used now the local specific enthalpy h = h{p,c) and with dh/dp = 1/p. 
Eq. ( [A16| ) means that if it is necessary to go beyond local thermodynamics, the non- 
dissipative (i.e., Hamiltonian or reactive) part of the stress tensor is modified according 
to 
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(A17) 



This can be written in a more intuitive form. The internal energy has the form 



8 



Hp, c) + fi'(Vc, c)] dr 



(A18) 



where the function g contains derivatives of c. Our example above had g{x, y) = K{y)x^/2. 
Now we define a generalized chemical potential, 



In terms of this chemical potential the body force is of the form 

— Vp + /icVc. 

This is the general form of the additional part of the non-dissipative part of the stress tensor 
used in the body of the paper. In the critical dynamics literature this is known as "Model 
H" in the Hohenberg-Halperin lexicon. Note that for an incompressible fluid one can 
rewrite this as 



where p' is an effective pressure chosen to guarantee V ■ v = 0. Intuitively gradients of the 
chemical potential generate body forces which need to be included in the non-dissipative 
part of the stress tensor. 

The equations presented here (when viscous dissipation is included) generalize the usual 
Navier-Stokes equation to situations which go beyond linear irreversible thermodynamics 
assumptions, but remain within the realm of classical fluid mechanics. 



/ic = 68 /6c. 



(A19) 



— Vp' — cV/ic 
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FIGURES 




Tt o oo vo 

O O O ON ON 



FIG. 1. Displacement of the drop's center of order parameter as a function of time for 
C = 1,R = 12, To = 1 and a = —0.004 in a square cell of side a = 200. Values for two sys- 
tem sizes are shown: o, a = 400; 0,0 = 200. The straight lines are linear fits to the data with 
slopes (dimensionless velocities) 0.0319 (a = 400) and 0.0316 (a = 200) respectively. 
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FIG. 3. Velocity of the center of order parameter as a function of the drop radius for 
C = l,ro = 0.5 in a square cell of side a = 200: o, a = -0.001; □, a = -0.004. Note that 
vcM is proportional to R for small values of R. Interestingly, linearity is lost at lower values of R 
at larger temperature gradients. 
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FIG. 4. Sequence of configurations showing the coalescence of two drops of radius i? = 12 in 
a cavity with a = b = 200. In this case C = l,a = —0.004 and tq = 1. The dimensionless times 
shown (from left to right and top to bottom) are t = 1305, 1555, 1605 and 2055. 
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FIG. 5. Spinodal decomposition of a binary fluid in a rectangular cefl of dimensions a = 800 
and b = 200 (vertical direction in the figures, corresponding to the y-direction in the fluid). The 
parameters used are C = 1, tq = 1 and a = —0.004. Note that the hotter end is at the bottom of 
each figure. Shown is the value of the order parameter in grey scale at three instants of time (from 
top to bottom): t = 6.25, 50 and 100. 
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